module cubed_sphere_array_mod

  use flogger
  use string
  use fiona
  use var_types_mod
  use cubed_sphere_objects
  use cubed_sphere_mesh_mod
  use cubed_sphere_domain_array_mod

  implicit none

  type cubed_sphere_array_type
    logical :: initialized = .false.
    type(var_stack_type) var_stack
    type(cubed_sphere_mesh_type), pointer :: mesh => null()
    type(cubed_sphere_domain_array_type), allocatable :: domain(:)
  contains
    procedure :: init            => cubed_sphere_array_init
    procedure :: clear           => cubed_sphere_array_clear
    procedure :: add_var         => cubed_sphere_array_add_var
    procedure :: var_idx         => cubed_sphere_array_var_idx
    procedure :: allocate_arrays => cubed_sphere_array_allocate_arrays
    generic   :: get_array       => cubed_sphere_array_get_array_2d_r8_1, &
                                    cubed_sphere_array_get_array_2d_r8_2, &
                                    cubed_sphere_array_get_array_3d_r8_1, &
                                    cubed_sphere_array_get_array_3d_r8_2, &
                                    cubed_sphere_array_get_array_4d_r8_1, &
                                    cubed_sphere_array_get_array_4d_r8_2, &
                                    cubed_sphere_array_get_array_5d_r8  , &
                                    cubed_sphere_array_get_array_2d_l_1 , &
                                    cubed_sphere_array_get_array_3d_l_1 , &
                                    cubed_sphere_array_get_array_2d_i4_1, &
                                    cubed_sphere_array_get_array_2d_i4_2, &
                                    cubed_sphere_array_get_array_3d_i4_1, &
                                    cubed_sphere_array_get_array_3d_i4_2
    procedure, private :: cubed_sphere_array_get_array_2d_r8_1
    procedure, private :: cubed_sphere_array_get_array_2d_r8_2
    procedure, private :: cubed_sphere_array_get_array_3d_r8_1
    procedure, private :: cubed_sphere_array_get_array_3d_r8_2
    procedure, private :: cubed_sphere_array_get_array_4d_r8_1
    procedure, private :: cubed_sphere_array_get_array_4d_r8_2
    procedure, private :: cubed_sphere_array_get_array_5d_r8
    procedure, private :: cubed_sphere_array_get_array_2d_l_1
    procedure, private :: cubed_sphere_array_get_array_3d_l_1
    procedure, private :: cubed_sphere_array_get_array_2d_i4_1
    procedure, private :: cubed_sphere_array_get_array_2d_i4_2
    procedure, private :: cubed_sphere_array_get_array_3d_i4_1
    procedure, private :: cubed_sphere_array_get_array_3d_i4_2
    procedure :: create_dataset  => cubed_sphere_array_create_dataset
    procedure :: append_dataset  => cubed_sphere_array_append_dataset
    procedure :: write           => cubed_sphere_array_write
    procedure :: write_only_init => cubed_sphere_array_write_only_init
    procedure :: create_init_dataset  => cubed_sphere_array_create_init_dataset
    procedure :: append_init_dataset  => cubed_sphere_array_append_init_dataset
    procedure :: write_init           => cubed_sphere_array_write_init
    final :: cubed_sphere_array_final
  end type cubed_sphere_array_type

contains

  subroutine cubed_sphere_array_init(this, mesh)

    class(cubed_sphere_array_type) , intent(inout), target :: this
    type(cubed_sphere_mesh_type)   , intent(in)   , target :: mesh

    integer idom

    call this%clear()

    this%mesh => mesh
    allocate(this%domain(mesh%ndom))
    do idom = 1, mesh%ndom
      if (this%mesh%domain(idom)%initialized) call this%domain(idom)%init(this%mesh%domain(idom))
    end do

    call this%var_stack%init()

    this%initialized = .true.

  end subroutine cubed_sphere_array_init

  subroutine cubed_sphere_array_clear(this)

    class(cubed_sphere_array_type), intent(inout) :: this

    integer idom

    call this%var_stack%clear()

    if (allocated(this%domain)) deallocate(this%domain)

    this%initialized = .false.

  end subroutine cubed_sphere_array_clear

  subroutine cubed_sphere_array_add_var(this, name, dtype, long_name, units, loc, tag, with_halo, fill_halo, output, only_2d, only_init)

    class(cubed_sphere_array_type), intent(inout) :: this
    character(*), intent(in) :: name
    character(*), intent(in), optional :: dtype
    character(*), intent(in), optional :: long_name
    character(*), intent(in), optional :: units
    character(*), intent(in), optional :: loc
    character(*), intent(in), optional :: tag
    character(*), intent(in), optional :: output
    character(*), intent(in), optional :: with_halo
    character(*), intent(in), optional :: fill_halo
    logical     , intent(in), optional :: only_2d
    logical     , intent(in), optional :: only_init

    call this%var_stack%append(name, dtype, long_name, units, loc, tag, with_halo, fill_halo, output, only_2d, only_init)

  end subroutine cubed_sphere_array_add_var

  integer function cubed_sphere_array_var_idx(this, name, loc) result(res)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in) :: name
    character(*), intent(in) :: loc

    do res = 1, this%var_stack%size
      if (this%var_stack%var_info(res)%name == name .and. this%var_stack%var_info(res)%loc == loc) return
    end do
    res = 0

  end function cubed_sphere_array_var_idx

  subroutine cubed_sphere_array_allocate_arrays(this)

    class(cubed_sphere_array_type), intent(inout) :: this

    integer idom
    logical need_update

    call this%var_stack%reorder()

    need_update = .true.
    do idom = 1, mesh%ndom
      if (this%domain(idom)%initialized) then
        call this%domain(idom)%allocate_arrays(this%var_stack, update_var_stack=need_update)
        need_update = .false.
      end if
    end do

  end subroutine cubed_sphere_array_allocate_arrays

  subroutine cubed_sphere_array_get_array_2d_r8_1(this, ptr, idom, ivar)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8), intent(out), pointer                     :: ptr(:,:)
    integer , intent(in )                              :: idom
    integer , intent(in )                              :: ivar

    if (.not. this%mesh%domain(idom)%initialized) then
      nullify(ptr)
      return
    end if
    call this%domain(idom)%get_array(ptr, this%var_stack%var_info(ivar))

  end subroutine cubed_sphere_array_get_array_2d_r8_1

  subroutine cubed_sphere_array_get_array_2d_r8_2(this, ptr, idom, var_name, loc)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8)    , intent(out), pointer                 :: ptr(:,:)
    integer     , intent(in )                          :: idom
    character(*), intent(in )                          :: var_name
    character(*), intent(in )                          :: loc

    integer ivar

    if (.not. this%mesh%domain(idom)%initialized) then
      nullify(ptr)
      return
    end if
    ivar = this%var_idx(var_name, loc)
    call this%domain(idom)%get_array(ptr, this%var_stack%var_info(ivar))

  end subroutine cubed_sphere_array_get_array_2d_r8_2

  subroutine cubed_sphere_array_get_array_3d_r8_1(this, ptr, idom, ivar)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8), intent(out), pointer                     :: ptr(:,:,:)
    integer , intent(in )                              :: idom
    integer , intent(in )                              :: ivar

    if (.not. this%mesh%domain(idom)%initialized) then
      nullify(ptr)
      return
    end if
    call this%domain(idom)%get_array(ptr, this%var_stack%var_info(ivar))

  end subroutine cubed_sphere_array_get_array_3d_r8_1

  subroutine cubed_sphere_array_get_array_3d_r8_2(this, ptr, idom, var_name, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8)    , intent(out), pointer                 :: ptr(:,:,:)
    integer     , intent(in )                          :: idom
    character(*), intent(in ), optional                :: var_name
    character(*), intent(in ), optional                :: loc
    character(*), intent(in ), optional                :: tag

    integer i, is, ie
    logical with_halo

    if (.not. this%mesh%domain(idom)%initialized) then
      nullify(ptr)
      return
    end if
    if (present(var_name) .and. present(loc)) then
      i = this%var_idx(var_name, loc)
      if (i == 0) then
        call log_error('Failed to get for ' // trim(var_name) // ' at ' // trim(loc) // '!')
      end if
      call this%domain(idom)%get_array(ptr, this%var_stack%var_info(i))
    else if (present(loc) .and. present(tag)) then
      call this%var_stack%get_idx(dtype_r8, loc, tag, is, ie, with_halo)
      call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)
      if (.not. associated(ptr)) then
        call log_error('Failed to get array for tag ' // trim(tag) // ' at ' // trim(loc) // '!')
      end if
    end if

  end subroutine cubed_sphere_array_get_array_3d_r8_2

  subroutine cubed_sphere_array_get_array_4d_r8_1(this, ptr, idom, ivar)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8), intent(out), pointer                     :: ptr(:,:,:,:)
    integer , intent(in )                              :: idom ! Domain index
    integer , intent(in )                              :: ivar ! Variable index

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%domain(idom)%get_array(ptr, this%var_stack%var_info(ivar))
    if (.not. associated(ptr)) then
      call log_error('Failed to get array for index ' // to_str(ivar) // '!')
    end if

  end subroutine cubed_sphere_array_get_array_4d_r8_1

  subroutine cubed_sphere_array_get_array_4d_r8_2(this, ptr, idom, var_name, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8)    , intent(out), pointer                 :: ptr(:,:,:,:)
    integer     , intent(in )                          :: idom      ! Domain index
    character(*), intent(in ), optional                :: var_name
    character(*), intent(in ), optional                :: loc       ! Variable location
    character(*), intent(in ), optional                :: tag       ! Variable tag

    integer i, is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    if (present(var_name) .and. present(loc)) then
      i = this%var_idx(var_name, loc)
      call this%domain(idom)%get_array(ptr, this%var_stack%var_info(i))
      if (.not. associated(ptr)) then
        call log_error('Failed to get for ' // trim(var_name) // ' at ' // trim(loc) // '!')
      end if
    else if (present(loc) .and. present(tag)) then
      call this%var_stack%get_idx(dtype_r8, loc, tag, is, ie, with_halo)
      call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)
      if (.not. associated(ptr)) then
        call log_error('Failed to get array for tag ' // trim(tag) // ' at ' // trim(loc) // '!')
      end if
    end if

  end subroutine cubed_sphere_array_get_array_4d_r8_2

  subroutine cubed_sphere_array_get_array_5d_r8(this, ptr, idom, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    real(r8)    , intent(out), pointer                 :: ptr(:,:,:,:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in )                          :: tag  ! Variable tag

    type(var_info_type), pointer :: var_info
    integer is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%var_stack%get_idx(dtype_r8, loc, tag, is, ie, with_halo)
    call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)

  end subroutine cubed_sphere_array_get_array_5d_r8

  subroutine cubed_sphere_array_get_array_2d_l_1(this, ptr, idom, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    logical     , intent(out), pointer                 :: ptr(:,:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in )                          :: tag  ! Variable tag

    type(var_info_type), pointer :: var_info
    integer is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%var_stack%get_idx(dtype_l, loc, tag, is, ie, with_halo)
    call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)

  end subroutine cubed_sphere_array_get_array_2d_l_1

  subroutine cubed_sphere_array_get_array_2d_i4_1(this, ptr, idom, loc, tag, var_name)

    class(cubed_sphere_array_type), intent(in), target :: this
    integer     , intent(out), pointer                 :: ptr(:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in ), optional                :: tag
    character(*), intent(in )                          :: var_name

    type(var_info_type), pointer :: var_info
    integer ivar
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    ivar = this%var_idx(var_name, loc)
    with_halo = this%var_stack%var_info(ivar)%with_halo
    call this%domain(idom)%get_array(ptr, loc, with_halo, this%var_stack%var_info(ivar))

  end subroutine cubed_sphere_array_get_array_2d_i4_1

  subroutine cubed_sphere_array_get_array_2d_i4_2(this, ptr, idom, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    integer     , intent(out), pointer                 :: ptr(:,:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in )                          :: tag  ! Variable tag

    type(var_info_type), pointer :: var_info
    integer is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%var_stack%get_idx(dtype_i4, loc, tag, is, ie, with_halo)
    call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)

  end subroutine cubed_sphere_array_get_array_2d_i4_2

  subroutine cubed_sphere_array_get_array_3d_l_1(this, ptr, idom, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    logical     , intent(out), pointer                 :: ptr(:,:,:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in )                          :: tag  ! Variable tag

    type(var_info_type), pointer :: var_info
    integer is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%var_stack%get_idx(dtype_l, loc, tag, is, ie, with_halo)
    call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)

  end subroutine cubed_sphere_array_get_array_3d_l_1

  subroutine cubed_sphere_array_get_array_3d_i4_1(this, ptr, idom, ivar)

    class(cubed_sphere_array_type), intent(in), target :: this
    integer, intent(out), pointer                      :: ptr(:,:,:)
    integer, intent(in )                               :: idom
    integer, intent(in )                               :: ivar

    if (.not. this%mesh%domain(idom)%initialized) then
      nullify(ptr)
      return
    end if
    call this%domain(idom)%get_array(ptr, this%var_stack%var_info(ivar))

  end subroutine cubed_sphere_array_get_array_3d_i4_1

  subroutine cubed_sphere_array_get_array_3d_i4_2(this, ptr, idom, loc, tag)

    class(cubed_sphere_array_type), intent(in), target :: this
    integer     , intent(out), pointer                 :: ptr(:,:,:,:)
    integer     , intent(in )                          :: idom ! Domain index
    character(*), intent(in )                          :: loc  ! Variable location
    character(*), intent(in )                          :: tag  ! Variable tag

    type(var_info_type), pointer :: var_info
    integer is, ie
    logical with_halo

    nullify(ptr)
    if (.not. this%mesh%domain(idom)%initialized) return
    call this%var_stack%get_idx(dtype_i4, loc, tag, is, ie, with_halo)
    call this%domain(idom)%get_array(ptr, loc, with_halo, is, ie)

  end subroutine cubed_sphere_array_get_array_3d_i4_2

  subroutine cubed_sphere_array_create_dataset(this, dataset_name, desc, file_path, file_prefix)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)                   :: dataset_name
    character(*), intent(in), optional         :: desc
    character(*), intent(in), optional         :: file_path
    character(*), intent(in), optional         :: file_prefix

    character(4) dims_eta(3), dims_xy(4), dims_z(5), dims_2d(3), dims_3d(4), hor_metric_dims(5), ver_metric_dims(6)
    integer ndx, ndy, nz, nv, npt, neq, nfq, nfqv, ncq, nvq

    call this%mesh%get_params(ndx=ndx, ndy=ndy, nz=nz, nv=nv, npt=npt, neq=neq, nfq=nfq, nfqv=nfqv, ncq=ncq, nvq=nvq)

    call fiona_init()

    dims_xy(1)     = 'np'
    dims_xy(2)     = 'i'
    dims_xy(3)     = 'j'
    dims_xy(4)     = 'dom'

    dims_z(1)      = 'np'
    dims_z(2)      = 'i'
    dims_z(3)      = 'j'
    dims_z(4)      = 'k'
    dims_z(5)      = 'dom'

    dims_eta(1)    = 'np'
    dims_eta(2)    = 'k'
    dims_eta(3)    = 'dom'

    dims_2d(1)     = 'i'
    dims_2d(2)     = 'j'
    dims_2d(3)     = 'dom'
    
    dims_3d(1)     = 'i'
    dims_3d(2)     = 'j'
    dims_3d(3)     = 'k'
    dims_3d(4)     = 'dom'
    
    hor_metric_dims(1) = 'rowh'
    hor_metric_dims(2) = 'colh'
    hor_metric_dims(3) = 'i'
    hor_metric_dims(4) = 'j'
    hor_metric_dims(5) = 'dom'
    
    ver_metric_dims(1) = 'rowv'
    ver_metric_dims(2) = 'colv'
    ver_metric_dims(3) = 'i'
    ver_metric_dims(4) = 'j'
    ver_metric_dims(5) = 'k'
    ver_metric_dims(6) = 'dom'

    if (.not. fiona_has_dataset(dataset_name)) then
      call fiona_create_dataset(dataset_name, desc=desc, file_path=file_path, file_prefix=file_prefix, mpi_comm=proc%comm)
      call fiona_add_dim(dataset_name, 'i'  , 'grid index in x-dir', '-', ndx, add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'j'  , 'grid index in y-dir', '-', ndy, add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'k'  , 'grid index in z-dir', '-', nz , add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'dom' , size=6, decomp=.true.)
      call fiona_add_dim(dataset_name, 'rowh', size=2)
      call fiona_add_dim(dataset_name, 'colh', size=2)
      call fiona_add_dim(dataset_name, 'rowv', size=2)
      call fiona_add_dim(dataset_name, 'colv', size=1)
      call fiona_add_dim(dataset_name, 'time', add_var=.true.)

      call fiona_add_dim(dataset_name, 'vtx', size=nv   )
      call fiona_add_dim(dataset_name, 'np' , size=npt  )
      call fiona_add_dim(dataset_name, 'eq' , size=nfq  )
      call fiona_add_dim(dataset_name, 'fq' , size=nfq  )
      call fiona_add_dim(dataset_name, 'fqv', size=nfqv )
      call fiona_add_dim(dataset_name, 'cq' , size=nvq  )
      call fiona_add_dim(dataset_name, 'vq' , size=nvq  )
      
      call fiona_add_var(dataset_name, 'x'     , 'x coordinate on each point'             , '-'            , dims_xy    , 'r8')
      call fiona_add_var(dataset_name, 'y'     , 'y coordinate on each point'             , '-'            , dims_xy    , 'r8')
      call fiona_add_var(dataset_name, 'z'     , 'z coordinate on each point'             , 'm'            , dims_z     , 'r8')
      call fiona_add_var(dataset_name, 'eta'   , 'eta coordinate on each point'           , '-'            , dims_eta   , 'r8')
      call fiona_add_var(dataset_name, 'xi'    , 'xi coordinate on each point'            , 'm'            , dims_eta   , 'r8')
      
      call fiona_add_var(dataset_name, 'lon', 'Longitude'                                 , 'degrees_east' , dims_2d    , 'r8')
      call fiona_add_var(dataset_name, 'lat', 'Latitude'                                  , 'degrees_north', dims_2d    , 'r8')
      
      call fiona_add_var(dataset_name,  'Ah', 'Horizonatal transformation matrix'         , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name, 'iAh', 'Horizonatal inversed transformation matrix', '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Gh', 'Horizonatal metric tensor'                 , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name, 'iGh', 'Horizonatal inversed metric tensor'        , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Jh', 'Horizonatal Jacobian (det(Gh))'            , '-'            , dims_2d        , 'r8')
      
      call fiona_add_var(dataset_name, 'iGv', 'Vertical inversed metric tensor'           , '-'            , ver_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Jv', 'Vertical Jacobian (det(Gv))'               , '-'            , dims_3d        , 'r8')

      call fiona_add_var(dataset_name, 'Jab', 'Jacobian (det(G))'                         , '-'            , dims_3d        , 'r8')
    end if

    call this%append_dataset(dataset_name)

  end subroutine cubed_sphere_array_create_dataset

  subroutine cubed_sphere_array_append_dataset(this, dataset_name)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)                   :: dataset_name

    character(4) dims_2d(4), dims_2d_eq(5), dims_2d_cq(5), dims_2d_fqv(5)
    character(4) dims_3d(5), dims_3d_eq(6), dims_3d_cq(6), dims_3d_fqv(6)
    character(4) dims_2d_init(3), dims_2d_eq_init(4), dims_2d_cq_init(4), dims_2d_fqv_init(4)
    character(4) dims_3d_init(4), dims_3d_eq_init(5), dims_3d_cq_init(5), dims_3d_fqv_init(5)
    character(30) var_name
    integer i

    dims_2d(1) = 'i'
    dims_2d(2) = 'j'
    dims_2d(3) = 'dom'
    dims_2d(4) = 'time'
    
    dims_3d(1) = 'i'
    dims_3d(2) = 'j'
    dims_3d(3) = 'k'
    dims_3d(4) = 'dom'
    dims_3d(5) = 'time'
    
    dims_2d_eq(1) = 'eq'    ; dims_2d_cq(1) = 'cq'   ; dims_2d_fqv(1) = 'fqv'  
    dims_2d_eq(2) = 'i'     ; dims_2d_cq(2) = 'i'    ; dims_2d_fqv(2) = 'i'    
    dims_2d_eq(3) = 'j'     ; dims_2d_cq(3) = 'j'    ; dims_2d_fqv(3) = 'j'    
    dims_2d_eq(4) = 'dom'   ; dims_2d_cq(4) = 'dom'  ; dims_2d_fqv(4) = 'dom'  
    dims_2d_eq(5) = 'time'  ; dims_2d_cq(5) = 'time' ; dims_2d_fqv(5) = 'time' 
    
    dims_3d_eq(1) = 'fq'    ; dims_3d_cq(1) = 'vq'   ; dims_3d_fqv(1) = 'fqv'  
    dims_3d_eq(2) = 'i'     ; dims_3d_cq(2) = 'i'    ; dims_3d_fqv(2) = 'i'   
    dims_3d_eq(3) = 'j'     ; dims_3d_cq(3) = 'j'    ; dims_3d_fqv(3) = 'j'   
    dims_3d_eq(4) = 'k'     ; dims_3d_cq(4) = 'k'    ; dims_3d_fqv(4) = 'k'   
    dims_3d_eq(5) = 'dom'   ; dims_3d_cq(5) = 'dom'  ; dims_3d_fqv(5) = 'dom' 
    dims_3d_eq(6) = 'time'  ; dims_3d_cq(6) = 'time' ; dims_3d_fqv(6) = 'time'

    dims_2d_init(1) = 'i'
    dims_2d_init(2) = 'j'
    dims_2d_init(3) = 'dom'
    
    dims_3d_init(1) = 'i'
    dims_3d_init(2) = 'j'
    dims_3d_init(3) = 'k'
    dims_3d_init(4) = 'dom'
    
    dims_2d_eq_init(1) = 'eq'    ; dims_2d_cq_init(1) = 'cq'   ; dims_2d_fqv_init(1) = 'fqv'  
    dims_2d_eq_init(2) = 'i'     ; dims_2d_cq_init(2) = 'i'    ; dims_2d_fqv_init(2) = 'i'    
    dims_2d_eq_init(3) = 'j'     ; dims_2d_cq_init(3) = 'j'    ; dims_2d_fqv_init(3) = 'j'    
    dims_2d_eq_init(4) = 'dom'   ; dims_2d_cq_init(4) = 'dom'  ; dims_2d_fqv_init(4) = 'dom'  
    
    dims_3d_eq_init(1) = 'fq'    ; dims_3d_cq_init(1) = 'vq'   ; dims_3d_fqv_init(1) = 'fqv'  
    dims_3d_eq_init(2) = 'i'     ; dims_3d_cq_init(2) = 'i'    ; dims_3d_fqv_init(2) = 'i'   
    dims_3d_eq_init(3) = 'j'     ; dims_3d_cq_init(3) = 'j'    ; dims_3d_fqv_init(3) = 'j'   
    dims_3d_eq_init(4) = 'k'     ; dims_3d_cq_init(4) = 'k'    ; dims_3d_fqv_init(4) = 'k'   
    dims_3d_eq_init(5) = 'dom'   ; dims_3d_cq_init(5) = 'dom'  ; dims_3d_fqv_init(5) = 'dom' 

    do i = 1, this%var_stack%size
      associate (info => this%var_stack%var_info(i))
      if (info%output) then
        if (info%loc == 'CA') then
          var_name = info%name
        else
          var_name = trim(info%name) // '_' // to_lower(info%loc)
        end if
        if (fiona_has_var(dataset_name, var_name)) cycle
        if (proc%is_root()) call log_notice('Define output for variable ' // trim(var_name) // '.', pid=proc%id)
        if (info%only_init) then
          if (info%only_2d) then
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_init, dtype_name(info%dtype))
            case ('CQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_cq_init, dtype_name(info%dtype))
            case ('LQ', 'RQ', 'BQ', 'TQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_eq_init, dtype_name(info%dtype))
            case ('UQ', 'DQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_fqv_init, dtype_name(info%dtype))
            end select
          else
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_init, dtype_name(info%dtype))
            case ('CQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_cq_init, dtype_name(info%dtype))
            case ('LQ', 'RQ', 'BQ', 'TQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_eq_init, dtype_name(info%dtype))
            case ('UQ', 'DQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_fqv_init, dtype_name(info%dtype))
            end select
          end if
        else
          if (info%only_2d) then
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d, dtype_name(info%dtype))
            case ('CQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_cq, dtype_name(info%dtype))
            case ('LQ', 'RQ', 'BQ', 'TQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_eq, dtype_name(info%dtype))
            case ('UQ', 'DQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_fqv, dtype_name(info%dtype))
            end select
          else
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d, dtype_name(info%dtype))
            case ('CQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_cq, dtype_name(info%dtype))
            case ('LQ', 'RQ', 'BQ', 'TQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_eq, dtype_name(info%dtype))
            case ('UQ', 'DQ')
              call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_fqv, dtype_name(info%dtype))
            end select
          end if
        end if
      end if
      end associate
    end do

  end subroutine cubed_sphere_array_append_dataset

  subroutine cubed_sphere_array_write_only_init(this, dataset_name, time_in_seconds, new_file, tag)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)           :: dataset_name
    real(8)     , intent(in), optional :: time_in_seconds
    logical     , intent(in), optional :: new_file
    character(*), intent(in), optional :: tag

    character(30) var_name
    integer , pointer :: ptr_3d_i4(:,:,:)
    real(r8), pointer :: ptr_2d_r8(:,:)
    real(r8), pointer :: ptr_3d_r8(:,:,:)
    real(r8), pointer :: ptr_4d_r8(:,:,:,:)
    integer ids, ide, nbx, jds, jde, nby, kds, kde, nz, dms, dme
    integer pc, pvs, pve, pes(6), pee(6), pqs, pqe
    integer nv, npt, neq, nfq, ncq, nfqv
    integer idom, ivar

    call this%mesh%get_params(ids=ids, ide=ide, nbx=nbx, &
                              jds=jds, jde=jde, nby=nby, &
                              kds=kds, kde=kde, nz=nz  , &
                              dms=dms, dme=dme,          &
                              pc=pc  , pvs=pvs, pve=pve, &
                              pqs=pqs, pqe=pqe,          &
                              pes=pes, pee=pee,          &
                              nv=nv, npt=npt,            &
                              neq=neq, ncq=ncq,          &
                              nfq=nfq, nfqv=nfqv)

    call fiona_start_output(dataset_name, time_in_seconds, new_file, tag)
    
    do idom = dms, dme
      if (this%mesh%domain(idom)%active) then
        do ivar = 1, this%var_stack%size
          associate (info => this%var_stack%var_info(ivar))
          if (info%output .and. info%only_init) then
            if (info%loc == 'CA') then
              var_name = info%name
            else
              var_name = trim(info%name) // '_' // to_lower(info%loc)
            end if
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              if (info%only_2d) then
                call this%get_array(ptr_2d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_2d_r8(ids:ide,jds:jde), start=[ids,jds,idom], count=[nbx,nby,1])
              else
                select case (info%dtype)
                case (dtype_r8)
                  call this%get_array(ptr_3d_r8, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_r8(ids:ide,jds:jde,kds:kde), start=[ids,jds,kds,idom], count=[nbx,nby,nz,1])
                case (dtype_i4)
                  call this%get_array(ptr_3d_i4, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_i4(ids:ide,jds:jde,kds:kde), start=[ids,jds,kds,idom], count=[nbx,nby,nz,1])
                end select
              end if
            case ('CQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[ncq,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[ncq,nbx,nby,nz,1])
              end if
            case ('LQ', 'RQ', 'BQ', 'TQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[nfq,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[nfq,nbx,nby,nz,1])
              end if
            case ('UQ', 'DQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[nfqv,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[nfqv,nbx,nby,nz,1])
              end if
            end select
          end if
          end associate
        end do
        
        associate (domain => this%domain(idom)%domain)
        call fiona_output(dataset_name, 'x'  , domain%x  (:,ids:ide,jds:jde        ), start=[1,ids,jds,    idom], count=[npt,nbx,nby,   1])
        call fiona_output(dataset_name, 'y'  , domain%y  (:,ids:ide,jds:jde        ), start=[1,ids,jds,    idom], count=[npt,nbx,nby,   1])
        call fiona_output(dataset_name, 'z'  , domain%z  (:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[npt,nbx,nby,nz,1])
        call fiona_output(dataset_name, 'eta', domain%eta(:,                kds:kde), start=[1,        kds,idom], count=[npt,        nz,1])
        call fiona_output(dataset_name, 'xi' , domain%xi (:,                kds:kde), start=[1,        kds,idom], count=[npt,        nz,1])
        
        call fiona_output(dataset_name, 'lon', domain%lon(    pc,ids:ide,jds:jde) * deg, start=[    ids,jds,idom], count=[    nbx,nby,1])
        call fiona_output(dataset_name, 'lat', domain%lat(    pc,ids:ide,jds:jde) * deg, start=[    ids,jds,idom], count=[    nbx,nby,1])
        call fiona_output(dataset_name, 'Ah' , domain% Ah(:,:,pc,ids:ide,jds:jde)      , start=[1,1,ids,jds,idom], count=[2,2,nbx,nby,1])
        call fiona_output(dataset_name, 'iAh', domain%iAh(:,:,pc,ids:ide,jds:jde)      , start=[1,1,ids,jds,idom], count=[2,2,nbx,nby,1])
        call fiona_output(dataset_name, 'Gh' , domain% Gh(:,:,pc,ids:ide,jds:jde)      , start=[1,1,ids,jds,idom], count=[2,2,nbx,nby,1])
        call fiona_output(dataset_name, 'iGh', domain%iGh(:,:,pc,ids:ide,jds:jde)      , start=[1,1,ids,jds,idom], count=[2,2,nbx,nby,1])
        call fiona_output(dataset_name, 'Jh' , domain% Jh(    pc,ids:ide,jds:jde)      , start=[    ids,jds,idom], count=[    nbx,nby,1])
        
        call fiona_output(dataset_name, 'iGv', domain%iGv(:,:,pc,ids:ide,jds:jde,kds:kde), start=[1,1,ids,jds,kds,idom], count=[2,1,nbx,nby,nz,1])
        call fiona_output(dataset_name, 'Jv' , domain% Jv(    pc,ids:ide,jds:jde,kds:kde), start=[    ids,jds,kds,idom], count=[    nbx,nby,nz,1])

        call fiona_output(dataset_name, 'Jab', domain%Jab(    pc,ids:ide,jds:jde,kds:kde), start=[    ids,jds,kds,idom], count=[    nbx,nby,nz,1])
        end associate
      endif
    end do
    call fiona_end_output(dataset_name)

  end subroutine cubed_sphere_array_write_only_init

  subroutine cubed_sphere_array_write(this, dataset_name, time_in_seconds, new_file, tag)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)           :: dataset_name
    real(8)     , intent(in), optional :: time_in_seconds
    logical     , intent(in), optional :: new_file
    character(*), intent(in), optional :: tag

    character(30) var_name
    integer , pointer :: ptr_3d_i4(:,:,:)
    real(r8), pointer :: ptr_2d_r8(:,:)
    real(r8), pointer :: ptr_3d_r8(:,:,:)
    real(r8), pointer :: ptr_4d_r8(:,:,:,:)
    integer ids, ide, nbx, jds, jde, nby, kds, kde, nz, dms, dme
    integer pc, pvs, pve, pes(6), pee(6), pqs, pqe
    integer nv, npt, neq, nfq, ncq, nfqv
    integer idom, ivar

    call this%mesh%get_params(ids=ids, ide=ide, nbx=nbx, &
                              jds=jds, jde=jde, nby=nby, &
                              kds=kds, kde=kde, nz=nz  , &
                              dms=dms, dme=dme,          &
                              pc=pc  , pvs=pvs, pve=pve, &
                              pqs=pqs, pqe=pqe,          &
                              pes=pes, pee=pee,          &
                              nv=nv, npt=npt,            &
                              neq=neq, ncq=ncq,          &
                              nfq=nfq, nfqv=nfqv)

    call fiona_start_output(dataset_name, time_in_seconds, new_file, tag)
    
    do idom = dms, dme
      if (this%mesh%domain(idom)%active) then
        do ivar = 1, this%var_stack%size
          associate (info => this%var_stack%var_info(ivar))
          if (info%output .and. .not. info%only_init) then
            if (info%loc == 'CA') then
              var_name = info%name
            else
              var_name = trim(info%name) // '_' // to_lower(info%loc)
            end if
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              if (info%only_2d) then
                call this%get_array(ptr_2d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_2d_r8(ids:ide,jds:jde), start=[ids,jds,idom], count=[nbx,nby,1])
              else
                select case (info%dtype)
                case (dtype_r8)
                  call this%get_array(ptr_3d_r8, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_r8(ids:ide,jds:jde,kds:kde), start=[ids,jds,kds,idom], count=[nbx,nby,nz,1])
                case (dtype_i4)
                  call this%get_array(ptr_3d_i4, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_i4(ids:ide,jds:jde,kds:kde), start=[ids,jds,kds,idom], count=[nbx,nby,nz,1])
                end select
              end if
            case ('CQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[ncq,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[ncq,nbx,nby,nz,1])
              end if
            case ('LQ', 'RQ', 'BQ', 'TQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[nfq,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[nfq,nbx,nby,nz,1])
              end if
            case ('UQ', 'DQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ids:ide,jds:jde), start=[1,ids,jds,idom], count=[nfqv,nbx,nby,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ids:ide,jds:jde,kds:kde), start=[1,ids,jds,kds,idom], count=[nfqv,nbx,nby,nz,1])
              end if
            end select
          end if
          end associate
        end do
      endif
    end do
    call fiona_end_output(dataset_name)

  end subroutine cubed_sphere_array_write

  subroutine cubed_sphere_array_create_init_dataset(this, dataset_name, desc, file_path, file_prefix)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)                   :: dataset_name
    character(*), intent(in), optional         :: desc
    character(*), intent(in), optional         :: file_path
    character(*), intent(in), optional         :: file_prefix

    character(4) dims_eta(3), dims_xy(4), dims_z(5), dims_2d(4), dims_3d(5), hor_metric_dims(6), ver_metric_dims(7), cs_dims(7)
    integer nx, ny, ndx, ndy, nz, nv, npt, neq, nfq, nfqv, ncq, nvq

    call this%mesh%get_params(nx=nx, ny=ny, ndx=ndx, ndy=ndy, nz=nz, nv=nv, npt=npt, neq=neq, nfq=nfq, nfqv=nfqv, ncq=ncq, nvq=nvq)

    call fiona_init()

    dims_xy(1)     = 'np'
    dims_xy(2)     = 'i'
    dims_xy(3)     = 'j'
    dims_xy(4)     = 'dom'

    dims_z(1)      = 'np'
    dims_z(2)      = 'i'
    dims_z(3)      = 'j'
    dims_z(4)      = 'k'
    dims_z(5)      = 'dom'

    dims_eta(1)    = 'np'
    dims_eta(2)    = 'k'
    dims_eta(3)    = 'dom'

    dims_2d(1)     = 'np'
    dims_2d(2)     = 'i'
    dims_2d(3)     = 'j'
    dims_2d(4)     = 'dom'
    
    dims_3d(1)     = 'np'
    dims_3d(2)     = 'i'
    dims_3d(3)     = 'j'
    dims_3d(4)     = 'k'
    dims_3d(5)     = 'dom'
    
    hor_metric_dims(1) = 'rowh'
    hor_metric_dims(2) = 'colh'
    hor_metric_dims(3) = 'np'
    hor_metric_dims(4) = 'i'
    hor_metric_dims(5) = 'j'
    hor_metric_dims(6) = 'dom'
    
    ver_metric_dims(1) = 'rowv'
    ver_metric_dims(2) = 'colv'
    ver_metric_dims(3) = 'np'
    ver_metric_dims(4) = 'i'
    ver_metric_dims(5) = 'j'
    ver_metric_dims(6) = 'k'
    ver_metric_dims(7) = 'dom'
    
    cs_dims(1) = 'cs1'
    cs_dims(2) = 'cs2'
    cs_dims(3) = 'cs3'
    cs_dims(4) = 'np'
    cs_dims(5) = 'id'
    cs_dims(6) = 'jd'
    cs_dims(7) = 'dom'

    if (.not. fiona_has_dataset(dataset_name)) then
      call fiona_create_dataset(dataset_name, desc=desc, file_path=file_path, file_prefix=file_prefix, mpi_comm=proc%comm)
      call fiona_add_dim(dataset_name, 'i'  , 'grid index in x-dir'       , '-', nx , add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'j'  , 'grid index in y-dir'       , '-', ny , add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'id' , 'domain grid index in x-dir', '-', ndx, add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'jd' , 'domain grid index in y-dir', '-', ndy, add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'k'  , 'grid index in z-dir'       , '-', nz , add_var=.false., decomp=.true.)
      call fiona_add_dim(dataset_name, 'dom' , size=6, decomp=.true.)
      call fiona_add_dim(dataset_name, 'rowh', size=2)
      call fiona_add_dim(dataset_name, 'colh', size=2)
      call fiona_add_dim(dataset_name, 'rowv', size=2)
      call fiona_add_dim(dataset_name, 'colv', size=1)
      call fiona_add_dim(dataset_name, 'cs1' , size=3)
      call fiona_add_dim(dataset_name, 'cs2' , size=3)
      call fiona_add_dim(dataset_name, 'cs3' , size=3)
      call fiona_add_dim(dataset_name, 'time', add_var=.true.)

      call fiona_add_dim(dataset_name, 'vtx', size=nv   )
      call fiona_add_dim(dataset_name, 'np' , size=npt  )
      call fiona_add_dim(dataset_name, 'eq' , size=nfq  )
      call fiona_add_dim(dataset_name, 'fq' , size=nfq  )
      call fiona_add_dim(dataset_name, 'fqv', size=nfqv )
      call fiona_add_dim(dataset_name, 'cq' , size=nvq  )
      call fiona_add_dim(dataset_name, 'vq' , size=nvq  )
      
      call fiona_add_var(dataset_name, 'x'     , 'x coordinate on each point'             , '-'            , dims_xy    , 'r8')
      call fiona_add_var(dataset_name, 'y'     , 'y coordinate on each point'             , '-'            , dims_xy    , 'r8')
      call fiona_add_var(dataset_name, 'z'     , 'z coordinate on each point'             , 'm'            , dims_z     , 'r8')
      call fiona_add_var(dataset_name, 'eta'   , 'eta coordinate on each point'           , '-'            , dims_eta   , 'r8')
      call fiona_add_var(dataset_name, 'xi'    , 'xi coordinate on each point'            , 'm'            , dims_eta   , 'r8')
      
      call fiona_add_var(dataset_name, 'lon', 'Longitude'                                 , 'degrees_east' , dims_2d    , 'r8')
      call fiona_add_var(dataset_name, 'lat', 'Latitude'                                  , 'degrees_north', dims_2d    , 'r8')
      
      call fiona_add_var(dataset_name,  'Ah', 'Horizonatal transformation matrix'         , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name, 'iAh', 'Horizonatal inversed transformation matrix', '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Gh', 'Horizonatal metric tensor'                 , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name, 'iGh', 'Horizonatal inversed metric tensor'        , '-'            , hor_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Jh', 'Horizonatal Jacobian (det(Gh))'            , '-'            , dims_2d        , 'r8')
      call fiona_add_var(dataset_name,  'CS', 'Christoffel symbol'                        , '-'            , cs_dims        , 'r8')
      
      call fiona_add_var(dataset_name, 'iGv', 'Vertical inversed metric tensor'           , '-'            , ver_metric_dims, 'r8')
      call fiona_add_var(dataset_name,  'Jv', 'Vertical Jacobian (det(Gv))'               , '-'            , dims_3d        , 'r8')

      call fiona_add_var(dataset_name, 'Jab', 'Jacobian (det(G))'                         , '-'            , dims_3d        , 'r8')
      
      !call fiona_add_var(dataset_name, 'lon_c' , 'Longitude of cell center'           , 'degrees_east' , [       'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_c' , 'Latitude of cell center'            , 'degrees_north', [       'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_v' , 'Longitude of cell vertices'         , 'degrees_east' , ['vtx', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_v' , 'Latitude of cell vertices'          , 'degrees_north', ['vtx', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_lq', 'Longitude of left edge'             , 'degrees_east' , ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_lq', 'Latitude of left edge'              , 'degrees_north', ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_rq', 'Longitude of right edge'            , 'degrees_east' , ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_rq', 'Latitude of right edge'             , 'degrees_north', ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_bq', 'Longitude of bottom edge'           , 'degrees_east' , ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_bq', 'Latitude of bottom edge'            , 'degrees_north', ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_tq', 'Longitude of top edge'              , 'degrees_east' , ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_tq', 'Latitude of top edge'               , 'degrees_north', ['eq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_dq', 'Longitude of down edge'             , 'degrees_east' , ['fqv', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_dq', 'Latitude of down edge'              , 'degrees_north', ['fqv', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_uq', 'Longitude of up edge'               , 'degrees_east' , ['fqv', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_uq', 'Latitude of up edge'                , 'degrees_north', ['fqv', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lon_cq', 'Longitude of cell quadrature points', 'degrees_east' , ['cq ', 'i  ', 'j  ', 'dom'], 'r8')
      !call fiona_add_var(dataset_name, 'lat_cq', 'Latitude of cell quadrature points' , 'degrees_north', ['cq ', 'i  ', 'j  ', 'dom'], 'r8')
    end if

    call this%append_dataset(dataset_name)

  end subroutine cubed_sphere_array_create_init_dataset

  subroutine cubed_sphere_array_append_init_dataset(this, dataset_name)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)                   :: dataset_name

    character(4) dims_2d(3), dims_2d_eq(5), dims_2d_cq(5), dims_2d_fqv(5)
    character(4) dims_3d(4), dims_3d_eq(6), dims_3d_cq(6), dims_3d_fqv(6)
    character(30) var_name
    integer i

    dims_2d(1) = 'i'
    dims_2d(2) = 'j'
    dims_2d(3) = 'dom'
    
    dims_3d(1) = 'i'
    dims_3d(2) = 'j'
    dims_3d(3) = 'k'
    dims_3d(4) = 'dom'
    
    dims_2d_eq(1) = 'eq'    ; dims_2d_cq(1) = 'cq'   ; dims_2d_fqv(1) = 'fqv'  
    dims_2d_eq(2) = 'i'     ; dims_2d_cq(2) = 'i'    ; dims_2d_fqv(2) = 'i'    
    dims_2d_eq(3) = 'j'     ; dims_2d_cq(3) = 'j'    ; dims_2d_fqv(3) = 'j'    
    dims_2d_eq(4) = 'dom'   ; dims_2d_cq(4) = 'dom'  ; dims_2d_fqv(4) = 'dom'  
    
    dims_3d_eq(1) = 'fq'    ; dims_3d_cq(1) = 'vq'   ; dims_3d_fqv(1) = 'fqv'  
    dims_3d_eq(2) = 'i'     ; dims_3d_cq(2) = 'i'    ; dims_3d_fqv(2) = 'i'   
    dims_3d_eq(3) = 'j'     ; dims_3d_cq(3) = 'j'    ; dims_3d_fqv(3) = 'j'   
    dims_3d_eq(4) = 'k'     ; dims_3d_cq(4) = 'k'    ; dims_3d_fqv(4) = 'k'   
    dims_3d_eq(5) = 'dom'   ; dims_3d_cq(5) = 'dom'  ; dims_3d_fqv(5) = 'dom' 

    do i = 1, this%var_stack%size
      associate (info => this%var_stack%var_info(i))
      if (info%output) then
        if (info%loc == 'CA') then
          var_name = info%name
        else
          var_name = trim(info%name) // '_' // to_lower(info%loc)
        end if
        if (fiona_has_var(dataset_name, var_name)) cycle
        if (proc%is_root()) call log_notice('Define output for variable ' // trim(var_name) // '.', pid=proc%id)
        if (info%only_2d) then
          select case (info%loc)
          case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d, dtype_name(info%dtype))
          case ('CQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_cq, dtype_name(info%dtype))
          case ('LQ', 'RQ', 'BQ', 'TQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_eq, dtype_name(info%dtype))
          case ('UQ', 'DQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_2d_fqv, dtype_name(info%dtype))
          end select
        else
          select case (info%loc)
          case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d, dtype_name(info%dtype))
          case ('CQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_cq, dtype_name(info%dtype))
          case ('LQ', 'RQ', 'BQ', 'TQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_eq, dtype_name(info%dtype))
          case ('UQ', 'DQ')
            call fiona_add_var(dataset_name, var_name, info%long_name, info%units, dims_3d_fqv, dtype_name(info%dtype))
          end select
        end if
      end if
      end associate
    end do

  end subroutine cubed_sphere_array_append_init_dataset

  subroutine cubed_sphere_array_write_init(this, dataset_name, time_in_seconds, new_file, tag)

    class(cubed_sphere_array_type), intent(in) :: this
    character(*), intent(in)           :: dataset_name
    real(8)     , intent(in), optional :: time_in_seconds
    logical     , intent(in), optional :: new_file
    character(*), intent(in), optional :: tag

    character(30) var_name
    integer , pointer :: ptr_3d_i4(:,:,:)
    real(r8), pointer :: ptr_2d_r8(:,:)
    real(r8), pointer :: ptr_3d_r8(:,:,:)
    real(r8), pointer :: ptr_4d_r8(:,:,:,:)
    integer ims, ime, jms, jme, ids, ide, jds, jde, nx, ny, ndx, ndy, nbx, nby, kds, kde, nz, dms, dme
    integer pc, pvs, pve, pes(6), pee(6), pqs, pqe
    integer nv, npt, neq, nfq, ncq, nfqv
    integer idom, ivar

    call this%mesh%get_params(ims=ims, ime=ime, ndx=ndx,          &
                              jms=jms, jme=jme, ndy=ndy,          &
                              ids=ids, ide=ide, nbx=nbx,          &
                              jds=jds, jde=jde, nby=nby,          &
                              kds=kds, kde=kde, nz=nz,            &
                              dms=dms, dme=dme,                   &
                              pc=pc  , pvs=pvs, pve=pve,          &
                              pqs=pqs, pqe=pqe,                   &
                              pes=pes, pee=pee,                   &
                              nv=nv, npt=npt,                     &
                              neq=neq, ncq=ncq,                   &
                              nfq=nfq, nfqv=nfqv)
    nx = ime - ims + 1
    ny = jme - jms + 1
    
    call fiona_start_output(dataset_name, time_in_seconds, new_file, tag)
    
    do idom = dms, dme
      if (this%mesh%domain(idom)%active) then
        do ivar = 1, this%var_stack%size
          associate (info => this%var_stack%var_info(ivar))
          if (info%output) then
            if (info%loc == 'CA') then
              var_name = info%name
            else
              var_name = trim(info%name) // '_' // to_lower(info%loc)
            end if
            select case (info%loc)
            case ('C', 'CA', 'L', 'R', 'B', 'T', 'U', 'D')
              if (info%only_2d) then
                call this%get_array(ptr_2d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_2d_r8(ims:ime,jms:jme), start=[ids,jds,idom], count=[nx,ny,1])
              else
                select case (info%dtype)
                case (dtype_r8)
                  call this%get_array(ptr_3d_r8, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_r8(ims:ime,jms:jme,kds:kde), start=[ids,jds,kds,idom], count=[nx,ny,nz,1])
                case (dtype_i4)
                  call this%get_array(ptr_3d_i4, idom, ivar)
                  call fiona_output(dataset_name, var_name, ptr_3d_i4(ims:ime,jms:jme,kds:kde), start=[ids,jds,kds,idom], count=[nx,ny,nz,1])
                end select
              end if
            case ('CQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ims:ime,jms:jme), start=[1,ids,jds,idom], count=[ncq,nx,ny,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ims:ime,jms:jme,kds:kde), start=[1,ids,jds,kds,idom], count=[ncq,nx,ny,nz,1])
              end if
            case ('LQ', 'RQ', 'BQ', 'TQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ims:ime,jms:jme), start=[1,ids,jds,idom], count=[nfq,nx,ny,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ims:ime,jms:jme,kds:kde), start=[1,ids,jds,kds,idom], count=[nfq,nx,ny,nz,1])
              end if
            case ('UQ', 'DQ')
              if (info%only_2d) then
                call this%get_array(ptr_3d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_3d_r8(:,ims:ime,jms:jme), start=[1,ids,jds,idom], count=[nfqv,nx,ny,1])
              else
                call this%get_array(ptr_4d_r8, idom, ivar)
                call fiona_output(dataset_name, var_name, ptr_4d_r8(:,ims:ime,jms:jme,kds:kde), start=[1,ids,jds,kds,idom], count=[nfqv,nx,ny,nz,1])
              end if
            end select
          end if
          end associate
        end do
        
        associate (domain => this%domain(idom)%domain)
        call fiona_output(dataset_name, 'x'  , domain%x  (:,ims:ime,jms:jme        ), start=[1,ids,jds,    idom], count=[npt,nx,ny,   1])
        call fiona_output(dataset_name, 'y'  , domain%y  (:,ims:ime,jms:jme        ), start=[1,ids,jds,    idom], count=[npt,nx,ny,   1])
        call fiona_output(dataset_name, 'z'  , domain%z  (:,ims:ime,jms:jme,kds:kde), start=[1,ids,jds,kds,idom], count=[npt,nx,ny,nz,1])
        call fiona_output(dataset_name, 'eta', domain%eta(:,                kds:kde), start=[1,        kds,idom], count=[npt,      nz,1])
        call fiona_output(dataset_name, 'xi' , domain%xi (:,                kds:kde), start=[1,        kds,idom], count=[npt,      nz,1])
        
        call fiona_output(dataset_name, 'lon', domain%lon(      :,ims:ime,jms:jme) * deg, start=[      1,ids,jds,idom], count=[      npt,nx,ny,1])
        call fiona_output(dataset_name, 'lat', domain%lat(      :,ims:ime,jms:jme) * deg, start=[      1,ids,jds,idom], count=[      npt,nx,ny,1])
        call fiona_output(dataset_name, 'Ah' , domain% Ah(:,:,  :,ims:ime,jms:jme)      , start=[1,1,  1,ids,jds,idom], count=[2,2,  npt,nx,ny,1])
        call fiona_output(dataset_name, 'iAh', domain%iAh(:,:,  :,ims:ime,jms:jme)      , start=[1,1,  1,ids,jds,idom], count=[2,2,  npt,nx,ny,1])
        call fiona_output(dataset_name, 'Gh' , domain% Gh(:,:,  :,ims:ime,jms:jme)      , start=[1,1,  1,ids,jds,idom], count=[2,2,  npt,nx,ny,1])
        call fiona_output(dataset_name, 'iGh', domain%iGh(:,:,  :,ims:ime,jms:jme)      , start=[1,1,  1,ids,jds,idom], count=[2,2,  npt,nx,ny,1])
        call fiona_output(dataset_name, 'Jh' , domain% Jh(      :,ims:ime,jms:jme)      , start=[      1,ids,jds,idom], count=[      npt,nx,ny,1])
        call fiona_output(dataset_name, 'CS' , domain% CS(:,:,:,:,ids:ide,jds:jde)      , start=[1,1,1,1,ids,jds,idom], count=[3,3,3,ncq,ndx,ndy,1])
        
        call fiona_output(dataset_name, 'iGv', domain%iGv(:,:,:,ims:ime,jms:jme,kds:kde), start=[1,1,1,ids,jds,kds,idom], count=[2,1,npt,nx,ny,nz,1])
        call fiona_output(dataset_name, 'Jv' , domain% Jv(    :,ims:ime,jms:jme,kds:kde), start=[    1,ids,jds,kds,idom], count=[    npt,nx,ny,nz,1])

        call fiona_output(dataset_name, 'Jab', domain%Jab(    :,ims:ime,jms:jme,kds:kde), start=[    1,ids,jds,kds,idom], count=[    npt,nx,ny,nz,1])

        !call fiona_output(dataset_name, 'lon_c' , domain%lon(pc           ,ids:ide,jds:jde) * deg, start=[  ids,jds,idom], count=[     ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_c' , domain%lat(pc           ,ids:ide,jds:jde) * deg, start=[  ids,jds,idom], count=[     ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_v' , domain%lon(pvs   :pve   ,ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[ nv ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_v' , domain%lat(pvs   :pve   ,ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[ nv ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_lq', domain%lon(pes(1):pee(1),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_lq', domain%lat(pes(1):pee(1),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_rq', domain%lon(pes(2):pee(2),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_rq', domain%lat(pes(2):pee(2),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_bq', domain%lon(pes(3):pee(3),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_bq', domain%lat(pes(3):pee(3),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_tq', domain%lon(pes(4):pee(4),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_tq', domain%lat(pes(4):pee(4),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_dq', domain%lon(pes(5):pee(5),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfqv,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_dq', domain%lat(pes(5):pee(5),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfqv,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_uq', domain%lon(pes(6):pee(6),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfqv,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_uq', domain%lat(pes(6):pee(6),ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[nfqv,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lon_cq', domain%lon(pqs   :pqe   ,ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[ncq ,ndx,ndy,1])
        !call fiona_output(dataset_name, 'lat_cq', domain%lat(pqs   :pqe   ,ids:ide,jds:jde) * deg, start=[1,ids,jds,idom], count=[ncq ,ndx,ndy,1])
        end associate
      endif
    end do
    call fiona_end_output(dataset_name)

  end subroutine cubed_sphere_array_write_init

  subroutine cubed_sphere_array_final(this)

    type(cubed_sphere_array_type), intent(inout) :: this

    call this%clear()

  end subroutine cubed_sphere_array_final

end module cubed_sphere_array_mod
